Obligatorio de Análisis Predictivo de Series Temporales¶

Conteo Vehícular y aproximación Poisson¶

Estudiantes:

Felipe Bastarrica (158687)

Emiliano Espíndola (187231)

Descripción¶

Se trabaja con datos reales obtenidos por la Intendencia de Montevideo, a través del Centro de Gestión de Movilidad, mediante sus contadores de autos ubicados en diferentes puntos de la ciudad. Estos datos vienen en intervalos acumulados de a 5 minutos y son captados por separado por cada carril que corre en un mismo sentido.

Los datos relevados haciendo uso de estos dispositivos, son de dominio público y la Intendencia los disponibiliza a través del sitio catálogo de datos que pertenece al estado.

El set de datos tiene muchas medidas con una frecuencia muy alta, muy útil y motivante su utilización.

Se pretende ajustar un modelo Poisson de media dependiente del tiempo para la cantidad de autos que transitan por una calle a lo largo de un día, obteniendo una distribución de probabilidad que nos permite simular estos eventos.

Formato de los datos originales¶

Originalmente los datos se presentan por archivos que acumulan períodos grandes de tiempo (meses, semestres, etc), variando según el año de relevamiento.

Cada uno de estos archivos posee, además del volumen de autos acumulados cada 5 minutos, el dato acumulado durante la última hora y otros parámetros que se muestran en la tabla:

| Dato | Tipo | Descripción | |-------|---------|-------------| | id_detector | Entero | Indica el id de la cámara que está monitoerando un determinado carril donde realiza el análisis de imagen para detectar la presencia de un vehículo | | id_carril | Entero | Es el número del carril que se está monitoreando. Va de 1 en adelante, siendo 1 el carril de más a la izquierda | | fecha | AAAA-MM-DD | Fecha en que fue tomada la muestra | | hora | hh:mm:ss | Hora en que fue tomada la muestra | | dsc_avenida | Texto | Nombre de la vía en la que se mide el tránsito | | dsc_int_anterior | Texto | Nombre de la vía que forma el cruce desde donde vienen los vehículos | | dsc_int_siguiente | Texto | Nombre de la vía que forma el cruce donde está el medidor, en general el mismo se encuentra un poco antes de esta vía. El sentido de circulación será desde el curce con dsc_int_anterior hacia el curce con dsc_int_siguiente. | | latitud | Float | latitud de donde está el medidor | | longitud | Float | longitud de donde está el medidor | | volumen | Enterno | Cantidad de vehículos detectados en el carril en los últimos 5 minutos | | volumen_hora | Enterno | Cantidad de vehículos detectados en el carril en la última hora |

Pre procesamiento¶

Para el trabajo obligatorio se realizaron varias transformaciones previas sobre los datos:

  • Se tomó en cuenta la zona de Sarmiento y Av. José Requena y García con dirección a la rambla. Por lo que se tuvo que filtrar los archivos de varios meses, tomando solo los datos referidos a este punto.

  • Se agregan los datos de los 3 carriles allí existentes, contando únicamente los autos que pasaron en total.

  • Se tomó en cuenta la fecha y hora para lograr armar una serie de tiempo

  • Se completaron los datos faltantes con la media de la serie con las que se contaba originalmente.

Instalaciones y bibliotecas utilizadas¶

In [ ]:
# Instalaciones
#install.packages("data.table")
#install.packages("caret")
In [ ]:
#Bibliotecas
library(astsa)
library(forecast)
library(data.table) 

#Ajuste de tamaño de gráficas
options(repr.plot.width=15, repr.plot.height=8)

Datos¶

In [16]:
# Ejemplo de los datos crudos descargados de la web
df_ejemplo <- read.csv("./dataframe_ejemplo.csv", )
head(df_ejemplo)
A data.frame: 6 × 11
id_detectorid_carrilfechahoradsc_avenidadsc_int_anteriordsc_int_siguientelatitudlongitudvolumenvolumen_hora
<int><int><chr><chr><chr><chr><chr><dbl><dbl><int><int>
0132017-01-0100:00:00SarmientoRequenaRambla-34.91537-56.168360 0
1132017-01-0100:05:00SarmientoRequenaRambla-34.91537-56.168360 0
2132017-01-0100:10:00SarmientoRequenaRambla-34.91537-56.16836672
3132017-01-0100:15:00SarmientoRequenaRambla-34.91537-56.16836224
4132017-01-0100:20:00SarmientoRequenaRambla-34.91537-56.16836336
5132017-01-0100:25:00SarmientoRequenaRambla-34.91537-56.16836672

Se elige la avenida Sarmiento para esta tarea, donde el mismo cuenta con 3 detectores en total, cubriendo cada uno de sus carriles. Luego de un procesamiento, se llega a un dataframe con una columna de fechas date y otra con volumen. Este procesamiento se hace fuera de la Notebook con el script main.py, donde se saca la infomación de cada csv, en este caso, todo 2017.

In [17]:
# Datos crudos obtenidos cada 5 minutos de cada sensor de carril 
df_raw = read.csv("./2017_2018_raw.csv", sep =";")
head(df_raw) 
tail(df_raw)
A data.frame: 6 × 2
datevolumen
<chr><int>
11/1/17 00:000
21/1/17 00:000
31/1/17 00:000
41/1/17 00:050
51/1/17 00:050
61/1/17 00:050
A data.frame: 6 × 2
datevolumen
<chr><int>
30783731/12/17 23:500
30783831/12/17 23:503
30783931/12/17 23:500
30784031/12/17 23:550
30784131/12/17 23:553
30784231/12/17 23:550

Definición de frecuencias¶

In [18]:
freq_5min = 365 * 24 * 60 / 5
freq_5min
freq_15min = 365 * 24 * 60 / 15
freq_15min
freq_sem_15min = 365 * 24 * 60 / ( 3 * 5 * (365 / 7) )
freq_sem_15min
105120
35040
672

Primer estudio de los datos¶

Se convierte a una serie de tiempo los datos cargados, sumarizándolos cada 15 minutos.

In [19]:
# Generar una serie de tiempo con datos cada 5 minutos
df_ts_raw = ts(df_raw$volumen, frequency=freq_5min, start=c(2017,1))

# Agregar la serie anterior cada 15 minutos (sumando los volúmenes)
df_aggregate_raw = aggregate(df_ts_raw, nfrequency = freq_15min, FUN = sum)
df_aggregate_ts_raw = ts(df_aggregate_raw, frequency = freq_15min*3, start=c(2017,1))

# Plot de la TS
tsplot(df_aggregate_ts_raw, main="Serie original sumarizada cada 15 minutos")

En esta serie se puede ver que no alcanza el 2018, porque faltan datos en varios momentos del año, para corregir esto, se generó una serie anual con la frecuencia deseada y se hico un merge con los datos originales.

In [20]:
# Copia de dataframe original
originaldf <- df_raw

# Creación de nueva columna para new date
originaldf$minAsPOSIX <- as.POSIXct(originaldf$date, format="%d/%m/%y %H:%M",tz="GMT")

# Generación del vector para cada 5 minutos
ndays <- 365 # días a generar
minAsNumeric <-seq(0, 60 * 60 * 24 * ndays, by = 60 * 5) # de 0 a 365 días con un paso de 5 minutos

# Conversión de secuencia a posix (formato calendario)
minAsPOSIX <- as.POSIXct(minAsNumeric, origin = "2017-01-01", tz = "GMT")

# Columna con datos calendario para 1 año 
base_df = data.frame(minAsPOSIX)

# Merge de datos calendario completos con los datos de volumen original
df_merged <- merge(base_df, originaldf, all.x = TRUE, by="minAsPOSIX")

# Agregación de datos de volumen por los carriles (los que tengan misma fecha y hora los suma)
newdf = aggregate(df_merged$volumen, by = list(Date=df_merged$minAsPOSIX), FUN=sum)
colnames(newdf) <- c('date', 'volumen')
In [21]:
ts = ts(newdf$volumen, frequency=freq_5min, start=c(2017,1))
tsplot(ts, main="Serie original con NAs donde faltan datos")

Ahora se obtiene un nuevo dataframe donde en la primer columna están los datos *date* de fecha completos en todo el tiempo deseado sin huecos, con su correspondiente sumarización en la columna *volumen*.

En caso de faltar datos, quedan NA, luego serán sustituidos por la media de la serie a trabajar.

Se analiza en períodos más cortos para visualizar la forma, reconocer mejor los días, las semanas y los fines de semana.

NOTA: El 1/1/2017 fue domingo.

In [22]:
# Muestras por día (1 muestra cada 5 minutos x 2)
samples_per_day = (60 / 5) * 24
samples_per_week = samples_per_day * 7

# Días a tomar para análisis más pequeño
days_to_crop = 21
initial_week = 35

# Intervalo
ini = initial_week * samples_per_week
fin = ini + samples_per_day * days_to_crop
st = 2017 + ini/length(ts)

# Serie recortada
df_crop = ts[ini : fin]

# Promedio de serie cortada y sustituir
mean = mean(df_crop, na.rm=TRUE)
print(paste("Media:", mean))
df_crop[is.na(df_crop)] = mean
df_ts_crop = ts(df_crop, start = st, freq = freq_5min)
[1] "Media: 29.6276912885061"
In [23]:
# Plot de la serie cortada
tsplot(df_ts_crop, main="Serie con marcas semanales")

# Plot barras en las semanas
lp = 2
picos = rep(c( 0, rep(150, lp/2), 0 ,rep(NA, samples_per_week - lp - 4), 0 , rep(150, lp/2), 0 ), days_to_crop/7)
ts_picos = ts(picos, start = st, freq = freq_5min)
lines(ts_picos, lwd = 3, col = 4)

Ajustes¶

Serie a utilizar¶

Visto que la serie difiere mucho con las estaciones del año, se decide utilizar el último 30% de la serie (setiembre a diciembre aproximadamente).

In [24]:
# Intervalo
ini = initial_week * samples_per_week
fin = length(ts)
st = 2017 + ini/length(ts)

# Serie recortada
df_crop = ts[ini : fin]

# Media de serie cortada y sustituir los NAs
mean = mean(df_crop, na.rm=TRUE)
print(paste("Media:", mean))
df_crop[is.na(df_crop)] = mean
ts_crop = ts(df_crop, start = st, freq = freq_5min)
[1] "Media: 30.2471385934919"
In [25]:
# Plot de la serie cortada
tsplot(ts_crop, main="Serie original cortada con marcas semanales")

# Plot barras en las semanas
lp = 2
picos = rep(c( 0, rep(150, lp/2), 0 ,rep(NA, samples_per_week - lp - 4), 0 , rep(150, lp/2), 0 ), 52 - initial_week)
ts_picos = ts(picos, start = st, freq = freq_5min)
lines(ts_picos, lwd = 3, col = 4)
In [26]:
acf(ts_crop)
In [27]:
pacf(ts_crop)

Del ACF se puede observar que hay una correlación alta con valores pasados, pero el PACF nos dice cuánto para atrás deberíamos de ir para ajustar la serie.

Periodicidad¶

In [28]:
# Estudio de la periodicidad 
samples_per_day = 288  # 24 * 60 / 5

s = ts_crop
s = s-mean(s) # Se quita la media
n = length(s)
I = abs(fft(s))^2 # FFT y módulo cuadrado
I = I[1:floor(n/2)] # Recorto el vector a las frecuencias observables
P = (4/n^2)*I # Escalado del periodograma
f = (0:(n/2-1))/n*(freq_5min) # Vector de frecuencias para hacer el gráfico (se multiplica por la frecuencia de la serie, para normalizar

# Ploteo del periodograma
plot(f, P, type="b", xlab="Frecuencia", main="Periodograma escalado", col=4, lwd=1, pch=19)
In [29]:
plot(f[0:1000], P[0:1000], type="b", xlab="Frecuencia", main="Periodograma escalado", col=4, lwd=1, pch=19)

Del periodograma se pueden extraer las frecuencias principales. Buscamos los máximos del vector P para obtener el valor de las frecuencias.

In [30]:
samples_per_day = 1

# Primera frecuencia de mayor importancia
f1 = f[which.max(P)]*samples_per_day
print(paste("f1:",f1))
print(paste("which.max(P1):",which.max(P)))

# Segunda frecuencia de mayor importancia
P[which.max(P)] = 0
f2 = f[which.max(P)]*samples_per_day
print(paste("f2:",f2))
print(paste("which.max(P2):",which.max(P)))

# Tercera frecuencia de mayor importancia
P[which.max(P)] = 0
f3 = f[which.max(P)]*samples_per_day
print(paste("f3:",f3))
print(paste("which.max(P3):",which.max(P)))

# Cuarta frecuencia de mayor importancia
P[which.max(P)] = 0
f4 = f[which.max(P)]*samples_per_day
print(paste("f4:",f4))
print(paste("which.max(P4):",which.max(P)))

# Quinta frecuencia de mayor importancia
P[which.max(P)] = 0
f5 = f[which.max(P)]*samples_per_day
print(paste("f5:",f5))
print(paste("which.max(P5):",which.max(P)))
[1] "f1: 364.978878537122"
[1] "which.max(P1): 121"
[1] "f2: 729.957757074243"
[1] "which.max(P2): 241"
[1] "f3: 1094.93663561137"
[1] "which.max(P3): 361"
[1] "f4: 416.684219663214"
[1] "which.max(P4): 138"
[1] "f5: 313.273537411029"
[1] "which.max(P5): 104"

Ajuste por perioricidad¶

In [115]:
t = time(ts_crop)
t_start_crop = t[1]
#f = 365

# Creación de fit con las 2 frecuencias mas significativas halladas
fit_df_ts_crop = lm(ts_crop ~ cos(2*pi*f1*t)+sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*f2*t)+cos(2*pi*f3*t)+sin(2*pi*f3*t))
#fit_df_ts_crop = lm(ts_crop ~ cos(2*pi*f*t)+sin(2*pi*f*t)+cos(4*pi*f*t)+sin(4*pi*f*t)+cos(6*pi*f*t)+sin(6*pi*f*t)) # en caso de mantener f fija

summary(fit_df_ts_crop)
Call:
lm(formula = ts_crop ~ cos(2 * pi * f1 * t) + sin(2 * pi * f1 * 
    t) + cos(2 * pi * f2 * t) + sin(2 * pi * f2 * t) + cos(2 * 
    pi * f3 * t) + sin(2 * pi * f3 * t))

Residuals:
    Min      1Q  Median      3Q     Max 
-62.319  -7.241  -0.529   7.595  76.616 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           30.24714    0.08896  340.02   <2e-16 ***
cos(2 * pi * f1 * t)  20.01699    0.12580  159.11   <2e-16 ***
sin(2 * pi * f1 * t)  -8.80723    0.12580  -70.01   <2e-16 ***
cos(2 * pi * f2 * t) -10.81332    0.12580  -85.95   <2e-16 ***
sin(2 * pi * f2 * t)  -2.34551    0.12580  -18.64   <2e-16 ***
cos(2 * pi * f3 * t)  -3.85440    0.12580  -30.64   <2e-16 ***
sin(2 * pi * f3 * t)   7.18327    0.12580   57.10   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.54 on 34555 degrees of freedom
Multiple R-squared:  0.5495,	Adjusted R-squared:  0.5494 
F-statistic:  7025 on 6 and 34555 DF,  p-value: < 2.2e-16

De el ajuste se puede saber que todas las componentes son significativas.

In [80]:
adjustment_fit_df_ts_crop = ts(fitted(fit_df_ts_crop), start=t_start_crop, freq=freq_5min)
tsplot(ts_crop, col=4, lwd=2, main = "Serie original vs ajuste lineal")
lines(adjustment_fit_df_ts_crop, lwd=2)

# Residuals del ajuste
residuals_fit = residuals(fit_df_ts_crop)
In [81]:
samples_per_day = 288
# Ajuste para las primeras 2 semanas del dataframe recortado
tsplot(head(ts_crop, samples_per_day*days_to_crop), col=4, lwd=2, main = "Principio de serie original vs ajuste lineal aumentado")
lines(head(adjustment_fit_df_ts_crop, samples_per_day*days_to_crop), lwd=2)
In [82]:
# Ajuste para las últimas 2 semanas del dataframe
tsplot(tail(ts_crop, samples_per_day*days_to_crop), col=4, lwd=2, main = "Final de serie original vs ajuste lineal aumentado")
lines(tail(adjustment_fit_df_ts_crop, samples_per_day*days_to_crop), lwd=2)
In [83]:
acf(residuals_fit)
In [84]:
pacf(residuals_fit)
In [85]:
checkresiduals(residuals_fit)
Warning message in modeldf.default(object):
“Could not find appropriate degrees of freedom for this model.”

Sigue habiendo una componente fuerte y 2 menores a pesar de haber ajustado la serie.

Auto Arima¶

In [86]:
# Auto ARIMA con Seasonal
autoarima = auto.arima(ts_crop, seasonal = TRUE)
ajuste_autoarima = ts(fitted(autoarima), start=t_start_crop, freq=freq_5min)
summary(autoarima)
residuals_autoarima = residuals(autoarima)
Series: ts_crop 
ARIMA(5,0,0) with non-zero mean 

Coefficients:
         ar1     ar2     ar3     ar4     ar5     mean
      0.4533  0.2501  0.1449  0.0760  0.0514  30.2471
s.e.  0.0054  0.0059  0.0060  0.0059  0.0054   1.6497

sigma^2 = 56.14:  log likelihood = -118643.6
AIC=237301.3   AICc=237301.3   BIC=237360.5

Training set error measures:
                        ME     RMSE      MAE  MPE MAPE MASE          ACF1
Training set -0.0002061926 7.491732 5.370428 -Inf  Inf  NaN -0.0001746736
In [87]:
acf(residuals_autoarima, na.action = na.pass)
In [88]:
pacf(residuals_autoarima, na.action = na.pass)

Con el Auto ARIMA se ve como el ACF y el PACF quedan limpios, sin encontrar correlación con valores pasados luego de aplicar un arima autoregresivo de orden 5, media móvil 0 e integrador 0.

Poisson¶

Se quiere armar una lista con vectores que representan el conteo de autos para determinado momento.

In [89]:
# Lista de vectores 
new_daily_vector = list()

# Muestras por día
samples_per_day_5min = 288

# Frecuencia a contemplar para el muestreo
days_to_sum = 7

# Muestras totales por días usados
samples_per_sum = samples_per_day_5min * days_to_sum

# Cuántos días se quiere muestrear
data_forward = 12

for (i in 1:data_forward){
    if (i != 1){
        new_daily_vector[[i]] = ts_crop[(i-1)*(samples_per_sum+1)]
    }else{
        new_daily_vector[[i]] = ts_crop[1]
    }
    for (j in 1:samples_per_sum){
        new_daily_vector[[i]] = append(new_daily_vector[[i]], round(ts_crop[(i-1)*samples_per_sum + j]),0) # Redondeo sino optim da Warning  
    } 
    new_daily_vector[[i]] = rev(new_daily_vector[[i]])
}
In [90]:
# Plot vectores diarios en misma grafica
ts1 = ts(new_daily_vector[[1]], start=0, freq=samples_per_day_5min)
tsplot(ts1, ylim = c(-1,130), main = "Superposición de días a modelar")

# Plot de los vectores hallados
for (i in 2:data_forward){
    ts = ts(new_daily_vector[[i]], start=0, freq=samples_per_day_5min)
    lines(ts, col = i)
}
In [91]:
# Transponer la matríz para tener en las filas las muestras para el mismo momento para la frecuencia determinada.
transposed_list = do.call(Map,c(f=c,new_daily_vector))

La distribución de Poisson se usa para modelar el número de eventos que ocurren en un proceso de Poisson. Sea $X∼P(λ)$, esto es, una variable aleatoria con distribución de Poisson donde el número medio de eventos que ocurren en un determinado intervalo con esta función de probabilidad:

$$P(X=x)= \frac{e^{-\lambda} \lambda^x}{x!}$$

Esto se obtiene con la función *dpois*. A esta se le ingresa la serie a ajustar, la función lambda a utilizar, y si se devuelve la probabilidad como su logaritmo. Al ser eventos independientes, las probabilidades se multiplican. Como se obtiene el logaritmo, éstas se terminan sumando y haciéndose negativas para hallar el mínimo ya que se quiere maximizar la verosimilitud.

In [92]:
# Lista para guardar los parámetros de la salida optim
est_vector_par = list()

# Lista para guardar los valores medios iterados
mean_vector_par = list()

for (i in 1:samples_per_sum){
    f <- function(lambda) {
        -sum(dpois(unlist(transposed_list[i]),lambda = lambda, log = TRUE), na.rm = TRUE)
    }
    est = optim(par = transposed_list[i], fn = f, method = "Brent", lower = 0, upper = 150)
    est_vector_par = append(est_vector_par, est$par)
    mean_vector_par = append(mean_vector_par, mean(as.numeric(unlist(transposed_list[i]))))
}

est_vector_par_ts = ts(est_vector_par, frequency = samples_per_day_5min, start=c(0,1))
mean_vector_par_ts = ts(mean_vector_par, frequency = samples_per_day_5min, start=c(0,1))
In [93]:
# Plot de vectores hallados
tsplot(est_vector_par_ts, col = 4, main = "Resultado de optimizador para cada punto")
lines(mean_vector_par_ts, col = 2)

El resultado de la optimización es la media de la cantidad de autos en un mismo tiempo para cada ventana semanal. Lo cual tiene sentido, porque el lambda es eso, la media de ocurrencias del evento.

In [94]:
# Ploteo de Lambda encontrado frente a los datos analizados
tsplot(est_vector_par_ts, col = 4, lwd = 5, ylim = c(-1,130), main = "Superposición vs evaluación para cada punto")

for (i in 2:data_forward){
    ts = ts(new_daily_vector[[i]], start=0, freq=samples_per_day_5min)
    lines(ts)
}
lines(est_vector_par_ts, col = 4, lwd = 2)

Ahora se quiere encontrar una función sinusoidal que pueda estimar la serie.

Se busca estimar, a través de la optimización simulando una distribución de Poisson, la función $\lambda(t)$ que "mejor" ajuste la media de cada evento para obtener dicha serie.

Se quieren encontrar los parámetros:

$$Params = \{a, b, c, d, e\}$$

Para la función:

$$\lambda(t) = a + b\ cos(2\pi f t) + c\ sin(2\pi f t) + d\ cos(4\pi f t) + e\ sin(4\pi f t)$$

Con $f = 365$ en este caso para no agregar complejidad a la función.

Definimos el vector de parámetros iniciales aproximados a los antes encontrados por el modelado lineal *lm* para darle menos trabajo al optimizador.

In [95]:
init.par = c(30, -20, -8, 4, -7, 3, 7)#,0)

x = round(ts_crop)
t = time(ts_crop)

f <- function(params) {
    lambdat = params[1] + params[2]*cos(2*pi*365*t) + params[3]*sin(2*pi*365*t) + params[4]*cos(4*pi*365*t) + params[5]*sin(4*pi*365*t)+ params[6]*cos(6*pi*365*t) + params[7]*sin(6*pi*365*t)
    #lambdat = params[1] + params[2]*cos(2*pi*365*(t+params[8])) + params[3]*sin(2*pi*365*(t+params[8])) + params[4]*cos(4*pi*365*(t+params[8])) + params[5]*sin(4*pi*365*(t+params[8]))+ params[6]*cos(6*pi*365*(t+params[8])) + params[7]*sin(6*pi*365*(t+params[8]))
    lambdat[which(lambdat < 0)] = 0 
    -sum(dpois(x, lambda = lambdat, log = TRUE))
}

est = optim(par = init.par, fn = f, method="BFGS", hessian=TRUE, control=list(trace=1, maxit=100))
est
initial  value 229029.552541 
iter  10 value 221087.821826
final  value 221083.133073 
converged
$par
  1. 30.2478677317613
  2. -20.4103127778594
  3. -7.70806472502441
  4. 2.0845725560699
  5. -9.78840590245098
  6. 0.472423245451734
  7. 7.38738164049264
$value
221083.133073162
$counts
function
53
gradient
12
$convergence
0
$message
NULL
$hessian
A matrix: 7 × 7 of type dbl
2345.29123 917.55915 952.31453-296.30262 914.3812-726.83419-53.95972
917.559151024.49427 457.19058 95.36251 449.1774-319.25510139.54150
952.31453 457.190581320.79690-503.13711 822.1967-774.83971 22.95248
-296.30262 95.36251-503.137111001.54179-317.6491 644.85112207.24128
914.38121 449.17738 822.19666-317.649091343.7494-745.07324272.70803
-726.83419-319.25510-774.83971 644.85112-745.07321413.67222 26.17709
-53.95972 139.54150 22.95248 207.24128 272.7080 26.17709931.61898

Luego de algunas iteraciones, se encuentran los valores de los parámetros de ajuste. A simple viste se ve una diferencia leve en amplitudes halladas respecto a la LM. Se grafica la serie $Lambda$ encontrada contra la serie original.

In [96]:
ajuste_f <- function(t) {
    lambdat = est$par[1] + est$par[2]*cos(2*pi*365*t) + est$par[3]*sin(2*pi*365*t) + est$par[4]*cos(4*pi*365*t) + est$par[5]*sin(4*pi*365*t) + + est$par[6]*cos(6*pi*365*t) + est$par[7]*sin(6*pi*365*t)
    #lambdat = est$par[1] + est$par[2]*cos(2*pi*365*(t+est$par[8])) + est$par[3]*sin(2*pi*365*(t+est$par[8])) + est$par[4]*cos(4*pi*365*(t+est$par[8])) + est$par[5]*sin(4*pi*365*(t+est$par[8])) + + est$par[6]*cos(6*pi*365*(t+est$par[8])) + est$par[7]*sin(6*pi*365*(t+est$par[8]))
}

eval_poisson = eval(ajuste_f(t))
tsplot(ts_crop, col=4, lwd=2, main = "Serie original truncada vs evaluación de Lambda(t)")
lines(eval_poisson, col = 2)
In [97]:
tsplot(head(ts_crop,288*7), col=4, lwd=2, main = "Zoom de serie original truncada vs evaluación de Lambda(t)")
lines(eval_poisson, col = 2, lwd = 3)
In [98]:
# Residuales encontrados por la resta de la serie original contra la evacluacion de Lambda
residuals_poisson = ts_crop - eval_poisson

tsplot(ts(residuals_fit, start = 0, freq=freq_5min), col = 4, main = "Residuales de la LM vs Residuales en Evaluación de Lambda(t)")
lines(ts(residuals_poisson, start = 0, freq=freq_5min), col = 2)
In [99]:
tsplot(head(ts(residuals_fit, start = 0, freq=freq_5min), 288*7), col = 4, main = "Head de Residuales de la LM vs Residuales en Evaluación de Lambda(t)")
lines(head(ts(residuals_poisson, start = 0, freq=freq_5min), 288*7), col = 2)
In [100]:
acf(residuals_poisson)
In [101]:
pacf(residuals_poisson)
In [102]:
# Varianza para cada Residuals
print("Varianza de Residuals Poisson")
var(residuals_poisson)
print("Varianza de ajuste lineal")
var(residuals_fit)
[1] "Varianza de Residuals Poisson"
278.306724772753
[1] "Varianza de ajuste lineal"
273.392003974394
In [103]:
# Root Mean Square Error para cada Residuals
print("RMSE de Residuals Poisson")
sqrt(mean(residuals_poisson^2))
print("RMSE de ajuste lineal")
sqrt(mean(residuals_fit^2))
[1] "RMSE de Residuals Poisson"
16.6822861884647
[1] "RMSE de ajuste lineal"
16.5343307631091

Las Varianzas y RMSE dan similares, un poco menores en el LM.

Cuando se incluyó una fase dentro de las componentes senoidales, se obtuvo un parámetro en la estimación muy chico, que a la hora de evaluar la serie la Varianza dió 282, un poco peor que sin este desfazaje. Esta prueba se hizo para analizar si, al agregar esta componente no lineal, se lograban mejores resutados.

Se grafican a continuación ambos modelados para compararlos visualmente.

In [104]:
tsplot(ts(fitted(fit_df_ts_crop), start = t_start_crop, freq = freq_5min), main = "Gráfico de la LM y la Evacluación de Lambda")
lines(eval_poisson, col = 2)
In [105]:
one_week_ts_crop = head(ts_crop,samples_per_day_5min*days_to_sum)

tsplot(one_week_ts_crop, main = "Gráfico de la LM y la Evacluación de Lambda")
lines(head(ts(fitted(fit_df_ts_crop), start = t_start_crop, freq = freq_5min),288*7), col = 4, lwd = 4)
lines(head(eval_poisson,samples_per_day_5min*days_to_sum), col = 2, lwd = 4)
legend(x = min(time(one_week_ts_crop), time(one_week_ts_crop)), y = max(one_week_ts_crop,one_week_ts_crop), legend=c("Serie original", "Linear Model", "Evaluación Poisson"), fill = c(1,4,2))

Ambos modelados son muy similares, dado que los parámetros hallados también los son.

Se ajustan los residuales de la función Lambda a un autoarima.

In [106]:
# Auto ARIMA para los residuals de la evalcuación de Poisson
autoarima_poisson = auto.arima(residuals_poisson)
autoarima_poisson_ts = ts(fitted(autoarima_poisson), start=t_start_crop, freq=freq_5min)
summary(autoarima_poisson)
residuals_autoarima = residuals(autoarima_poisson)
Series: residuals_poisson 
ARIMA(5,0,2) with non-zero mean 

Coefficients:
         ar1      ar2     ar3     ar4     ar5      ma1     ma2     mean
      1.3683  -0.8478  0.2187  0.0788  0.1403  -0.9452  0.6921  -0.0007
s.e.  0.0260   0.0485  0.0185  0.0093  0.0115   0.0249  0.0408   0.7064

sigma^2 = 54.14:  log likelihood = -118018.3
AIC=236054.6   AICc=236054.6   BIC=236130.7

Training set error measures:
                       ME     RMSE      MAE      MPE     MAPE MASE        ACF1
Training set 5.849429e-05 7.357473 5.292776 -12.5111 243.4675  NaN -0.00541282
In [107]:
acf(residuals_autoarima)
In [108]:
pacf(residuals_autoarima)
In [109]:
checkresiduals(residuals_autoarima)
Warning message in modeldf.default(object):
“Could not find appropriate degrees of freedom for this model.”

Conclusiones¶

  • Se hace un preprocesamiento de los datos obtenidos para guardarlos sumando los sensores para la intersección seleccionada
  • Se estudia la serie, viendo la forma de la misma y limpiándola para extraer una serie temporal interesante de estudiar
  • Pudo encontrarse, aplicando diferentes técnicas, modelos que logren ajustar la serie estudiada